#ifndef _QGCM_MEXAUX_SRC_
#define _QGCM_MEXAUX_SRC_
#include "qgcmio.h"
#include "mat.h"
#include "matrix.h"
#include "mexAux.h"

// ..................................................................................................................... //
static void mxFetchDomainDataFromModel(model_class & mo, mxStruct & mxDomainData){

	const unsigned int numfields = mo.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue();
	mxDomainData.setField( "NumFields", mxContainer( numfields) );

	double ft_coeff = static_cast<double>( mo.get_global_param_by_name(gdb_str_FOURIER_TRANSFORM_COEFFICIENT).getfValue());
	mxDomainData.setField( "FCOEFF", mxContainer( ft_coeff ) );

	// * write the non zero field index values of A */
	const vector<rank3> & aindexmap = mo.getAfieldindices();
	if(aindexmap.size() == 0){
		fprintf(stderr, "Warning: A operator has no nonzero sets of field indices.\n");
		fprintf(stderr, "There is nonlinear data to write out.\n");
	}else{
		mxContainer mxAIndexMap( aindexmap );
		mxDomainData.setField( "AFieldIndices", mxContainer( mxAIndexMap ) );
	}

	const char * profNames [] = { 
		gdb_prof_str_K_DOMAIN, 
		gdb_prof_str_DK_DOMAIN, 
		gdb_prof_str_K_DOMAIN_CENTERED, 
		gdb_prof_str_NK_DOMAIN, 
		gdb_prof_str_FCOEFFS, 
		NULL 
	};
	for(int i = 0; profNames[i] != NULL; i++){
		const parameter & prof = mo.get_global_prof_by_name( profNames[i] );
		mxDomainData.setField( profNames[i], mxContainer(prof) );
	}
};
// ..................................................................................................................... //
bool mx_write_static_weights(const bool nonlinear, const static_container * Rsc, const static_container * Ssc, 
		const map_container<REAL, 6> * Amap,
		const map_container<REAL, 3> * Qmap, const map_container<REAL, 3> * QImap,
		MATWriter & mfile)
{
	// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	string description;

	if(nonlinear){
		// call function twice
		mfile.mxWriteVariable( mxStruct( *Rsc, "Nonzero elements (rank11 pos, value) of the R static coefficients."), "R_static_weights" );
	}
	// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	if(nonlinear){
		mxStruct sWeights( *Ssc, "Nonzero elements (rank11 pos, value) of the S static coefficients." );
		mfile.mxWriteVariable( sWeights, "S_static_weights");
	}
	// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	mxStruct amap(*Amap, "Nonzero elements (rank6 pos, value) of the A operator.");
	mfile.mxWriteVariable( amap, "A_operator" );
	// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	mxStruct qmap(*Qmap, "Nonzero elements (rank3 pos, value) of the Q operator."); ;
	mfile.mxWriteVariable( qmap, "Q_operator" );
	// . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
	mxStruct qimap(*QImap, "Nonzero elements (rank3 pos, value) of the QI operator.");;
	mfile.mxWriteVariable( qimap, "QI_operator" );
	return true;
};
// ---------------------------------------------------------------------------------------
#ifdef DEBUG
void dump_static_weights(const rank11 * RSindexmap, const qgcomplex * Rs, const unsigned int map_size){
	FILE * fout = getrout();
	for(unsigned int u = 0; u < map_size; ++u){
		const rank11 & r11 = RSindexmap[u];
		fprintf(fout, "|");
		for(size_t x = 0; x < 11; x++){
			fprintf(fout, " %4i", static_cast<int>(r11[x]));
		}
		fprintf(fout, "|");
		fprintf(fout, " = ( %15g, %15g )\n", Rs[u].real(), Rs[u].imag());
	}
};
// ---------------------------------------------------------------------------------------
#endif
// ---------------------------------------------------------------------------------------
// Function below writes the static weight data to a designated MAT file.
// ---------------------------------------------------------------------------------------
bool mx_write_static_coeffs
(
 model_class & mo, const char * fname, bool overwrite,
 const static_container * Rsc, const static_container * Ssc,
 const map_container<REAL, 6> * Amap,
 const map_container<REAL, 3> * Qmap,
 const map_container<REAL, 3> * QImap
 ){
	bool nonlinear = mo.get_global_param_by_name(gdb_str_NONLINEAR).getiValue() != 0;
	MATWriter matfile( fname );
	if(!matfile.isValid()){
		fprintf(stderr, "Error: The static coefficients output MATLAB file '%s' already exists.\n", 
				fname);
		fprintf(stderr, "Cannot overwrite existing static coefficients output file.\n");
		fprintf(stderr, "Please rename or delete file or use 'mode=overwrite' switch.\n");
		return false;
	}

	mxStruct mxDomainData;
	mxFetchDomainDataFromModel(mo, mxDomainData);
	matfile.mxWriteVariable( mxDomainData, "DomainData" );
	// ......................................
	if(! mx_write_static_weights( nonlinear, Rsc, Ssc, Amap, Qmap, QImap, matfile) ){
		fprintf(stderr, "Error writing static weights to file.\n"); 
		return false;
	}
	return true;
};
// ---------------------------------------------------------------------------------------
static bool mxCompareStructFields( const mxStruct & struct1, const mxStruct & struct2, const char * fieldsToCompare [] ){
	for(int i = 0; fieldsToCompare[i] != NULL; i++){
		mxContainer struct1Field = struct1.getField(fieldsToCompare[i]);
		if( !struct2.doesFieldExist(fieldsToCompare[i]) ){
			fprintf(stderr, "Cannot find necessary domain data field '%s'.", fieldsToCompare[i]);	
			return false;
		}	
		mxContainer struct2Field = struct2.getField(fieldsToCompare[i]);
		if( !(struct1Field == struct2Field ) ) return false;
	}
	return true;
}
// ---------------------------------------------------------------------------------------
static bool mxReadStaticDomainData(model_class & mo, const MATReader & mfile, const char * filename){
	// first check to make sure there is agreement between the global attribute FOURIER_TRANSFORM_COEFFICIENT
	// and that defined by the model here.

	// Pull the domain data from the mat file.
	mxStruct mxMatFileDomainData = mfile.mxGetStruct("DomainData" );
	if( mxMatFileDomainData.isEmpty() ) {
		fprintf(stderr, "Could not find domain data in static coefficient file.");
		return false;
	}

	// Pull the domain data from the current model instantiation:
	mxStruct mxModelDomainData;
	mxFetchDomainDataFromModel(mo, mxModelDomainData);

	const char * fieldsToCompare[] = { "FCOEFF", "K_DOMAIN", "DK_DOMAIN", "K_DOMAIN_CENTERED", "NK_DOMAIN", "FCOEFFS", NULL };
	// List above must be NULL terminated!
	if( ! mxCompareStructFields( mxModelDomainData, mxMatFileDomainData, fieldsToCompare )){
		fprintf(stderr, "Error: Mismatch between current model instantiation and static coefficients input file.\n");
		return false;
	}	

	// All the data we want here should be packed into the struct read into mxDomainData
	mxContainer mxFCoeffs = mxMatFileDomainData.getField("FCOEFF");
	if( mxFCoeffs.isEmpty() || !mxIsDouble(mxFCoeffs.mxGetArray()) ){
		fprintf(stderr, "Could not find domain data for the fourier transform cofficients.");
		return false;
	}
	if( ! ( mxFCoeffs ==  mo.get_global_param_by_name(gdb_str_FOURIER_TRANSFORM_COEFFICIENT) ) ){
		fprintf(stderr, "Error: Mismatch between current model instantiation and static coefficients input file '%s' in the parameter '%s'.\n", 
				filename, gdb_str_FOURIER_TRANSFORM_COEFFICIENT );
		fprintf(stderr, "You may be inadvertently mixing fourier transform conventions.\n");
		return false;
	}
	// need to read the Aindexmap and store it

	mxContainer mxAfieldIndexMap( mfile.mxGetVariable("AFieldIndices") );
	if( ! mxAfieldIndexMap.isEmpty() ){
		const mxArray * mxIndexList = mxAfieldIndexMap.mxGetArray();
		const mwSize * dims = mxGetDimensions( mxIndexList );
		const int size = dims[0];
		rassert( mxGetNumberOfDimensions( mxIndexList ) == 2 );	
		const short * mxIndexListData = static_cast<short*>(mxGetData( mxIndexList ));
		vector<rank3> Aindexmap( size, rank3('*','*','*'));
		for(int row=0; row < size; row++){
			rank3 indices;
			for(int col=0; col < 3; col++){
				indices[col] = mxIndexListData[ col*size+row ];
			}
			Aindexmap[row] = indices;
		}	
		mo.setAfieldindices( Aindexmap ); // assign the nonzero index map to the model for general purposes
	};
	return true;
};

// ---------------------------------------------------------------------------------------
template <class T, int WIDTH>
static bool mxReadStaticVar(const char * varname, map_container<T, WIDTH> & target, const MATReader & mfile){
	mxStruct mxStaticVar = mfile.mxGetStruct( varname );
	string descr;
	if( !mxStaticVar.unpack( target, descr ) ){
		fprintf(stderr, "Unable to load static coefficients variable '%s'.\n", varname);
		return false;
	}
	return true;
}

template
static bool mxReadStaticVar(const char * varname, map_container<REAL, 3> & target, const MATReader & mfile);


// ---------------------------------------------------------------------------------------
bool mx_read_static_data(model_class & mo, const char * fname, map_container<REAL,3> &  Qcoeffs, map_container<REAL,3> & QIcoeffs){
	MATReader mfile( fname );
	if( ! mfile.isValid() ) {
		fprintf(stderr, "Unable to open MAT file '%s'.", fname);
		return false;
	}
	if( ! mxReadStaticDomainData( mo, mfile, fname ) ){
		return false;
	}
	if( ! mxReadStaticVar( "Q_operator", Qcoeffs, mfile ) ){
		return false;
	}
	if( ! mxReadStaticVar( "QI_operator", QIcoeffs, mfile ) ){
		return false;
	}
	return true;
};
// ---------------------------------------------------------------------------------------
bool mx_read_static_data(model_class & mo, const char * fname,
		static_container & Rs, static_container & Ss, map_container<REAL,6> & Acoeffs,
		map_container<REAL,3> &  Qcoeffs, map_container<REAL,3> & QIcoeffs){
	MATReader mfile( fname );
	if( ! mfile.isValid() ) {
		fprintf(stderr, "Unable to open MAT file '%s'.", fname);
		return false;
	}
	if( ! mxReadStaticVar( "R_static_weights" , Rs, mfile ) ){
		return false;
	}   
	if( ! mxReadStaticVar( "S_static_weights" , Ss, mfile ) ){
		return false;
	} 
	if( ! mxReadStaticVar( "A_operator", Acoeffs, mfile ) ){
		return false;
	}
	if( ! mxReadStaticVar( "Q_operator", Qcoeffs, mfile ) ){
		return false;
	}
	if( ! mxReadStaticVar( "QI_operator", QIcoeffs, mfile ) ){
		return false;
	}
	return true;
};
// ---------------------------------------------------------------------------------------
// Machinery for packing/unpacking the model state into a MAT file
// ---------------------------------------------------------------------------------------

typedef const char * const MXMATLABEL;

class mxModelInterface {
	private:
		vector<string> errorList;
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		class mxSaveModelStateFailure : std::exception {
			const char * what() {
				return "mxSaveModelState::Failure";
			}
		};
	protected:
		const char * const mfileName;
		bool success;

		static MXMATLABEL LABEL_MODEL_STATE;
		static MXMATLABEL LABEL_MODEL_PARAMETERS;
		static MXMATLABEL LABEL_MODEL_PROFILES;
		static MXMATLABEL LABEL_MODEL_PSYSTEM;	
		static MXMATLABEL LABEL_MODEL_SSYSTEM;
		static MXMATLABEL LABEL_MODEL_SCRIPT_INFO;	
		static MXMATLABEL LABEL_MODEL_SCRIPT_NAME ;
		static MXMATLABEL LABEL_MODEL_SCRIPT_ID ; 
		static MXMATLABEL LABEL_MODEL_SCRIPT_TYPE ;
		static MXMATLABEL LABEL_MODEL_SCRIPT_TYPE_FUNC ;
		static MXMATLABEL LABEL_MODEL_SCRIPT_TYPE_SCRIPT ;
		static MXMATLABEL LABEL_MODEL_SCRIPT_CALLSIG ;
		static MXMATLABEL LABEL_MODEL_SCRIPT_NUMARGS ;
		static MXMATLABEL LABEL_MODEL_SCRIPT_RETURNTYPE ;
		static MXMATLABEL LABEL_MODEL_SCRIPT_BODY;
		static MXMATLABEL LABEL_MODEL_DIMENSIONS;
		static MXMATLABEL LABEL_MODEL_DIMENSIONS_NUMBER_TIME_STEPS;
		static MXMATLABEL LABEL_MODEL_DIMENSIONS_NUMBER_OF_FIELDS;
		static MXMATLABEL LABEL_MODEL_DIMENSIONS_NUMBER_OF_ANGULAR_MODES;
		static MXMATLABEL LABEL_MODEL_DIMENSIONS_NUMBER_FOURIER_MODES;
		static MXMATLABEL LABEL_MODEL_PSYSTEM_VARIABLE_HISTORY;
		static MXMATLABEL LABEL_MODEL_PSYSTEM_MEANFIELD_HISTORY;
		static MXMATLABEL LABEL_MODEL_PSYSTEM_FUNCTION_EXPR;
		static MXMATLABEL LABEL_MODEL_PSYSTEM_FUNCTION_LAST_VALUE;
		static MXMATLABEL LABEL_MODEL_PSYSTEM_FUNCTION_LIST;
		static MXMATLABEL LABEL_MODEL_PSYSTEM_FUNCTION_DERIVATIVE_OF;
		static MXMATLABEL LABEL_MODEL_PSYSTEM_FLUXFIELD_INDICES;
		static MXMATLABEL LABEL_MODEL_PSYSTEM_FLUXFIELD_LAST_VALUE;
		static MXMATLABEL LABEL_MODEL_PSYSTEM_FLUXFIELD_LIST;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_DISSIPOP_INDEX;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_DISSIPOP_EXPR;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_DISSIPOP_LIST;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_INDEX_NK;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_FIELD_INDEX1;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_FIELD_INDEX2;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_LIST;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_EXPR;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_LINEAROP_INDEX_NK;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_LINEAROP_FIELD_INDEX1;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_LINEAROP_FIELD_INDEX2;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_LINEAROP_EXPR;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_LINEAROP_LIST;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_TIMESTEP_DESCRIPTION;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_TIMESTEP_DATA;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_CURRENT_TIMESTEP;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_TIMESTEP_LABEL;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_DATA;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_REGSET_INDICES;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_REGSET_PRIMARY_DATA;
		static MXMATLABEL LABEL_MODEL_SSYSTEM_REGSET_INTERMEDIATE_DATA;

	protected:
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		void lAssert(bool b, const string & errorOnFailure){
			if( ! b ) {
				errorList.push_back(errorOnFailure);
				throw mxSaveModelStateFailure(); 
			}
		}
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		void dumpErrors() {
			if( ! errorList.empty() ){
				fprintf(stderr, "Errors occurred while attempting to save the model state to MAT file '%s'.\n",
						mfileName);
				for(size_t i = 0; i < errorList.size(); i++){
					fprintf(stderr, "%s\n", errorList[i].c_str() );
				}
			}
		}
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		static string sFormatError( const char * fmt, ... ){
			va_list args;
			va_start(args, fmt);
			string rval;
			stdvsprintf(rval, fmt, args);
			return rval;
		}
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
	public:
		mxModelInterface(const char * fname) : success(false), mfileName(fname) {};

		bool operator==( const bool b ) const {
			return this->success == b;
		}
		bool didExecutionSucceed() const {
			return this->success;
		}
};

MXMATLABEL mxModelInterface::LABEL_MODEL_STATE = "ModelState";
MXMATLABEL mxModelInterface::LABEL_MODEL_PARAMETERS = "Parameters";
MXMATLABEL mxModelInterface::LABEL_MODEL_PROFILES = "Profiles";
MXMATLABEL mxModelInterface::LABEL_MODEL_SCRIPT_INFO = "Scripts";	
MXMATLABEL mxModelInterface::LABEL_MODEL_PSYSTEM = "PSystem";	
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM = "SSystem";	
MXMATLABEL mxModelInterface::LABEL_MODEL_SCRIPT_NAME = "Name";
MXMATLABEL mxModelInterface::LABEL_MODEL_SCRIPT_ID = "ID"; 
MXMATLABEL mxModelInterface::LABEL_MODEL_SCRIPT_TYPE = "ScriptType";
MXMATLABEL mxModelInterface::LABEL_MODEL_SCRIPT_TYPE_FUNC = "Function";
MXMATLABEL mxModelInterface::LABEL_MODEL_SCRIPT_TYPE_SCRIPT = "Script";
MXMATLABEL mxModelInterface::LABEL_MODEL_SCRIPT_CALLSIG = "CallSignature";
MXMATLABEL mxModelInterface::LABEL_MODEL_SCRIPT_NUMARGS = "NumberOfArguments";
MXMATLABEL mxModelInterface::LABEL_MODEL_SCRIPT_RETURNTYPE = "ReturnType";
MXMATLABEL mxModelInterface::LABEL_MODEL_SCRIPT_BODY = "Body";
MXMATLABEL mxModelInterface::LABEL_MODEL_DIMENSIONS = "Dimensions";
MXMATLABEL mxModelInterface::LABEL_MODEL_DIMENSIONS_NUMBER_TIME_STEPS = "time";
MXMATLABEL mxModelInterface::LABEL_MODEL_DIMENSIONS_NUMBER_OF_FIELDS = "fields";
MXMATLABEL mxModelInterface::LABEL_MODEL_DIMENSIONS_NUMBER_OF_ANGULAR_MODES = "fourier-angular";
MXMATLABEL mxModelInterface::LABEL_MODEL_DIMENSIONS_NUMBER_FOURIER_MODES = "fourier";
MXMATLABEL mxModelInterface::LABEL_MODEL_PSYSTEM_VARIABLE_HISTORY = "Variables";
MXMATLABEL mxModelInterface::LABEL_MODEL_PSYSTEM_MEANFIELD_HISTORY = "Fields";
MXMATLABEL mxModelInterface::LABEL_MODEL_PSYSTEM_FUNCTION_LIST = "Functions";
MXMATLABEL mxModelInterface::LABEL_MODEL_PSYSTEM_FUNCTION_EXPR = "Expression";
MXMATLABEL mxModelInterface::LABEL_MODEL_PSYSTEM_FUNCTION_LAST_VALUE = "CurrentValue";
MXMATLABEL mxModelInterface::LABEL_MODEL_PSYSTEM_FUNCTION_DERIVATIVE_OF = "DerivativeOf";
MXMATLABEL mxModelInterface::LABEL_MODEL_PSYSTEM_FLUXFIELD_INDICES = "Indices";
MXMATLABEL mxModelInterface::LABEL_MODEL_PSYSTEM_FLUXFIELD_LAST_VALUE = "CurrentValue";
MXMATLABEL mxModelInterface::LABEL_MODEL_PSYSTEM_FLUXFIELD_LIST = "FluxFields";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_DISSIPOP_INDEX = "FieldIndex";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_DISSIPOP_EXPR = "Expression";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_DISSIPOP_LIST = "DissipationOperators";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_INDEX_NK = "AngularModeIndex";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_FIELD_INDEX1 = "FieldIndex1";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_FIELD_INDEX2 = "FieldIndex2";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_EXPR = "InhomogeneousOperator";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_LIST = "InhomogeneousOperators";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_LINEAROP_INDEX_NK = "AngularModeIndex";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_LINEAROP_FIELD_INDEX1 = "FieldIndex1";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_LINEAROP_FIELD_INDEX2 = "FieldIndex-2";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_LINEAROP_EXPR = "LinearOperator";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_LINEAROP_LIST = "LinearOperators";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_TIMESTEP_DESCRIPTION = "Description";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_TIMESTEP_DATA = "TimeStepWidths";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_CURRENT_TIMESTEP = "CurrentTimeStep";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_TIMESTEP_LABEL = "Timesteps";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_DATA = "SSystemData";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_REGSET_INDICES = "Indices";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_REGSET_PRIMARY_DATA = "PrimaryData";
MXMATLABEL mxModelInterface::LABEL_MODEL_SSYSTEM_REGSET_INTERMEDIATE_DATA = "IntermediateData";

// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
class mxSaveModelState : public mxModelInterface {
	private :
		const model_class & mo;
		MATWriter mout;
		mxStruct modelState;	
	private:
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		void checkMatFileStatus() {
			lAssert( mout.isValid(), 
					sFormatError("Error opening file '%s' for writing.\n", mfileName)	
			       );
		}
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		void errorCleanUp() {
			// if there is anything to clean up before destruction	
		};
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		void saveModelParameters() {
			mxStruct modelParams;
			size_t numParams = mo.get_num_global_params();
			for(size_t i = 0; i < numParams; i++){
				const parameter & param = mo.get_global_param_by_index( i );
				mxContainer mxParam( param );
				modelParams.setField( param.GetName(), mxParam );			
			}
			modelState.setField( LABEL_MODEL_PARAMETERS, modelParams );
		}
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		void saveModelProfiles() {
			mxStruct modelProfs;
			size_t numProfs = mo.get_num_global_profs();
			for(size_t i = 0; i < numProfs; i++){
				const parameter & prof = mo.get_global_prof_by_index( i );
				mxContainer mxProf( prof );
				modelProfs.setField( prof.GetName(), mxProf );			
			}
			modelState.setField( LABEL_MODEL_PROFILES, modelProfs );
		}	
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		void saveModelScriptInfo() {
			vector<const model_diagnostics_block*> func_list;
			mo.getdiaglist( func_list ); // retrieve pointers to diag blocks that are functions
			size_t flistsize = func_list.size();
			mxCellArrayOfStructs scriptInfoList( flistsize );
			for(size_t iter = 0; iter < flistsize; iter++){ // iterate through all the functions, etc, write data 
				mxStruct & scriptBlockInfo = scriptInfoList[iter];
				const model_diagnostics_block & func = *(func_list[iter]);
				scriptBlockInfo.setField( LABEL_MODEL_SCRIPT_NAME, mxContainer( func.getfuncname() ) );
				scriptBlockInfo.setField( LABEL_MODEL_SCRIPT_ID, mxContainer( static_cast<int>(iter) ) );
				if( func.isblockfunction() ){
					const int numargs = func.getfuncnumargs();
					const vector<char_type> cstacktypes = func.GetCallStackTypes(); 
					const char_type rtype = func.GetReturnType();
					string args;
					create_function_call_stack_string( cstacktypes, args );
					scriptBlockInfo.setField( LABEL_MODEL_SCRIPT_TYPE, mxContainer( LABEL_MODEL_SCRIPT_TYPE_FUNC ) );
					scriptBlockInfo.setField( LABEL_MODEL_SCRIPT_CALLSIG, mxContainer( args ) );
					scriptBlockInfo.setField( LABEL_MODEL_SCRIPT_NUMARGS, mxContainer( numargs ) );
					scriptBlockInfo.setField( LABEL_MODEL_SCRIPT_RETURNTYPE, mxContainer( get_type_id(rtype) ) );
				}else{
					scriptBlockInfo.setField( LABEL_MODEL_SCRIPT_TYPE, mxContainer(LABEL_MODEL_SCRIPT_TYPE_SCRIPT));

				}
				string body = func.getBlockAsText();
				scriptBlockInfo.setField( LABEL_MODEL_SCRIPT_BODY, mxContainer( body ) );
			}
			modelState.setField( LABEL_MODEL_SCRIPT_INFO, scriptInfoList );

		}
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		void saveModelDimensions() {
			mxStruct modelDims;

			// add the relevant dimensions
			int tindex = mo.getPsystem().get_tindex() + 1;
			modelDims.setField( LABEL_MODEL_DIMENSIONS_NUMBER_TIME_STEPS, tindex );

			// number of fields
			int numfields = mo.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue(); 
			modelDims.setField( LABEL_MODEL_DIMENSIONS_NUMBER_OF_FIELDS, numfields);

			// number of fourier angular modes
			const vector<int> & nkdomain = ( mo.get_global_prof_by_name(gdb_prof_str_NK_DOMAIN).getiData() );
			modelDims.setField( LABEL_MODEL_DIMENSIONS_NUMBER_OF_ANGULAR_MODES, nkdomain.size() );

			// number of kcentered values
			const vector<REAL> & kdomain_centered = (mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData());
			modelDims.setField( LABEL_MODEL_DIMENSIONS_NUMBER_FOURIER_MODES, kdomain_centered.size() );

			modelState.setField( LABEL_MODEL_DIMENSIONS, modelDims );
		}
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		void saveModelPSystem() {
			// start by saving all variables in a time array
			mxStruct psystemInfo;
			const psystem & psys = mo.getPsystem();
			const vector< vector< Value > > & variableHistory = psys.getvariablehistory();
			const vector< vector< Value > > & meanFieldHistory = psys.getmeanfieldhistory();

			const size_t numtimes = psys.getnumtsteps(); // the last one is not stored in the history
			mxStruct mxVarHistoryList;
			for(size_t iter = 0; iter < psys.getvariablesize(); iter++){
				// start by creating the name for the variable field
				const string & varname = psys.getvariablename(iter);
				mxVector<double> mxVarHistory( numtimes, 0 );
				for(size_t ti = 0; ti < numtimes - 1; ti++){
					mxVarHistory[ti] = static_cast<double>( variableHistory[ti][iter].GetFloat());
				}
				const Variable & variable = psys.getvariablebyindex( iter );				
				mxVarHistory[numtimes - 1] = static_cast<double>( variable.GetFloat() );
				mxVarHistoryList.setField( varname.c_str(), mxVarHistory );
			}
			psystemInfo.setField( LABEL_MODEL_PSYSTEM_VARIABLE_HISTORY, mxVarHistoryList );

			mxStruct mxMeanfieldHistoryList;
			for(size_t iter = 0; iter < psys.getmeanfieldsize(); iter++){
				const string & fieldname = psys.getmeanfieldname(iter);
				mxVector<double> mxFieldHistory( numtimes, 0 );

				for(size_t ti = 0; ti < numtimes - 1; ti++){
					mxFieldHistory[ti] = static_cast<double>( meanFieldHistory[ti][iter].GetFloat() ); 
				}
				const Variable & meanfield = psys.getmeanfield( iter );
				mxFieldHistory[numtimes - 1] = static_cast<double>( meanfield.GetFloat() );
				mxMeanfieldHistoryList.setField( fieldname.c_str(), mxFieldHistory );
			}
			psystemInfo.setField( LABEL_MODEL_PSYSTEM_MEANFIELD_HISTORY, mxMeanfieldHistoryList );
			// now save function information (this must include it's ParserX expression ...)
			// as well as the an index to the mean field it points to if the function is a derivative
			// tricky stuff
			mxStruct mxFunctionList;
			for(size_t iter = 0; iter < psys.getfunctionssize(); iter++){
				mxStruct mxFuncInfo;
				mxFuncInfo.setField( LABEL_MODEL_PSYSTEM_FUNCTION_EXPR, psys.getfunctionexpr(iter) );
				const Variable & value = psys.getfunctionvalue(iter);
				mxFuncInfo.setField( LABEL_MODEL_PSYSTEM_FUNCTION_LAST_VALUE, value.GetFloat() );
				int dindex = -1;
				for( size_t x = 0; dindex < 0 && x < psys.getmeanfieldsize(); x++){
					const Variable & ptr = psys.getmeanfield(iter);
					if( value == ptr ) dindex = x;
				}
				if(dindex > 0){
					const string & mean_field_name = psys.getmeanfieldname(dindex);
					mxFuncInfo.setField( LABEL_MODEL_PSYSTEM_FUNCTION_DERIVATIVE_OF, mean_field_name );
				}
				mxFunctionList.setField( psys.getfunctionname(iter).c_str(), mxFuncInfo );
			}

			psystemInfo.setField( LABEL_MODEL_PSYSTEM_FUNCTION_LIST, mxFunctionList );
			// finally save the flux fields
			mxStruct mxFluxfields;
			for(size_t iter = 0; iter < psys.getsmeanfieldsize(); iter++){
				const rank2 & r2 = psys.getsmeanfieldindices(iter);
				const string & name = psys.getsmeanfieldname(iter);
				const Variable & val = psys.getsmeanfield(iter);
				mxVector<double> mxIndices(2,0);
				for(size_t i = 0; i < 2; i++) mxIndices[i] = static_cast<double>(r2[i]);
				mxStruct mxFluxfield;
				mxFluxfield.setField( LABEL_MODEL_PSYSTEM_FLUXFIELD_INDICES, mxIndices);
				mxFluxfield.setField( LABEL_MODEL_PSYSTEM_FLUXFIELD_LAST_VALUE, val.GetFloat() );
				mxFluxfields.setField( name.c_str(), mxFluxfield );
			}
			psystemInfo.setField( LABEL_MODEL_PSYSTEM_FLUXFIELD_LIST, mxFluxfields );

			modelState.setField( LABEL_MODEL_PSYSTEM , psystemInfo ); 
		}
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		void saveModelSSystem() {
			mxStruct ssystemInfo;

			// save dissipation operators first
			const ssystem & ssys = mo.getSsystem();
			const vector< pair<INDEX, ParserX> > & dissipation_ops = ssys.get_dissipation_ops();
			const vector< pair<rank3, ParserX> > & inhomogeneous_ops = ssys.get_inhomogeneous_ops();
			const lopinfo & linear_ops = ssys.get_linear_op_info();
			const int num_d_ops = static_cast<int>(dissipation_ops.size());

			mxCellArrayOfStructs mxDissipOps( dissipation_ops.size() );
			for(size_t iter = 0; iter < dissipation_ops.size(); iter++){
				mxStruct & mxDissipOp = mxDissipOps[iter];
				const INDEX & index = dissipation_ops[iter].first;
				const ParserX & p = dissipation_ops[iter].second;
				mxDissipOp.setField( LABEL_MODEL_SSYSTEM_DISSIPOP_INDEX,  static_cast<double>( index ) );
				mxDissipOp.setField( LABEL_MODEL_SSYSTEM_DISSIPOP_EXPR,  p.GetExprAsString().c_str() );
			}
			ssystemInfo.setField( LABEL_MODEL_SSYSTEM_DISSIPOP_LIST, mxDissipOps );

			mxCellArrayOfStructs mxInhomogeneousOps( inhomogeneous_ops.size() );
			for(size_t iter = 0; iter < inhomogeneous_ops.size(); iter++){
				mxStruct & mxInhomogeneousOp = mxInhomogeneousOps[iter];
				const rank3 & indices = inhomogeneous_ops[iter].first;
				const ParserX & p = inhomogeneous_ops[iter].second;
				mxInhomogeneousOp.setField( LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_INDEX_NK, static_cast<double>(indices[0]) );
				mxInhomogeneousOp.setField( LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_FIELD_INDEX1, static_cast<double>(indices[1]) );
				mxInhomogeneousOp.setField( LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_FIELD_INDEX2, static_cast<double>(indices[2]) );
				mxInhomogeneousOp.setField( LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_EXPR, p.GetExprAsString().c_str() );
			}
			ssystemInfo.setField( LABEL_MODEL_SSYSTEM_INHOMOGENEOUSOP_LIST, mxInhomogeneousOps );

			mxCellArrayOfStructs mxLinearOps( linear_ops.lop_list.size() );
			for( size_t iter = 0; iter < linear_ops.lop_list.size(); iter++) {
				const lopinfo::element & lop = linear_ops.lop_list[iter];
				mxStruct & mxLinearOp = mxLinearOps[iter];
				mxLinearOp.setField( LABEL_MODEL_SSYSTEM_LINEAROP_INDEX_NK, static_cast<double>(lop.angular_mode) );
				mxLinearOp.setField( LABEL_MODEL_SSYSTEM_LINEAROP_FIELD_INDEX1, static_cast<double>(lop.field_indices[0]) );
				mxLinearOp.setField( LABEL_MODEL_SSYSTEM_LINEAROP_FIELD_INDEX2, static_cast<double>(lop.field_indices[1]) );
				mxLinearOp.setField( LABEL_MODEL_SSYSTEM_LINEAROP_EXPR, lop.parser->GetExprAsString().c_str() );
			}
			ssystemInfo.setField( LABEL_MODEL_SSYSTEM_LINEAROP_LIST, mxLinearOps );

			mxStruct mxTimeStepInfo;

			mxTimeStepInfo.setField( LABEL_MODEL_SSYSTEM_TIMESTEP_DESCRIPTION, "The nondimensional timestep size at each time step.");
			mxVector<double> mxTimeSteps( ssys.get_tindex(), 0 );
			for(size_t iter = 0; iter < ssys.get_tindex(); iter++) mxTimeSteps[iter] = ssys.gettimestep(iter);
			mxTimeStepInfo.setField( LABEL_MODEL_SSYSTEM_TIMESTEP_DATA, mxTimeSteps);
			mxTimeStepInfo.setField( LABEL_MODEL_SSYSTEM_CURRENT_TIMESTEP, static_cast<double>( mo.get_global_param_by_name(gdb_str_TIMESTEP).getfValue()));
			ssystemInfo.setField( LABEL_MODEL_SSYSTEM_TIMESTEP_LABEL, mxTimeStepInfo );

			mxStruct mxSsystemData;

			saveSsystemDataSet( mxSsystemData, "phi", ssys.get_phi_data() );
			saveSsystemDataSet( mxSsystemData, "G", ssys.get_G_data() );
			saveSsystemDataSet( mxSsystemData, "Fphi", ssys.get_F_phi_data() );
			saveSsystemDataSet( mxSsystemData, "FG", ssys.get_F_G_data() );
			saveSsystemDataSet( mxSsystemData, "R", ssys.get_R_data() );
			saveSsystemDataSet( mxSsystemData, "S", ssys.get_S_data() );

			if( ssys.is_model_using_regularization() ){
				saveSsystemDataSet( mxSsystemData, "Fphi_regularized", ssys.get_F_phi_regularized_data(), true );
				saveSsystemDataSet( mxSsystemData, "R_regularized", ssys.get_R_regularized_data(), true );
				saveSsystemDataSet( mxSsystemData, "S_regularized", ssys.get_S_regularized_data(), true );
			}	
			ssystemInfo.setField( LABEL_MODEL_SSYSTEM_DATA, mxSsystemData );
			modelState.setField( LABEL_MODEL_SSYSTEM , ssystemInfo ); 

		}
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		static bool shouldWriteComplexData( const qgcomplex & data ){
			// here we check if the complex data is not close enough to zero to be ignored
			static const REAL zero(0);
			const REAL* dataPtrs[2] = { &(data.real()), &(data.imag()) };
			for(size_t i = 0; i < 2; i++){
				const int * ptr =  (int*) dataPtrs[i];
				int intDiff = *ptr - *(int*)(&zero);
			        if(intDiff < -4 || intDiff > 4) return true;
			}
			return false;
		}
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		void saveSsystemDataSet( mxStruct & mxDest, const char * fieldName, const twotimevar & dataSet, bool equalTimeOnly = false ) {
			const int maxTime = dataSet.getmaxtime();
			/* Data sets are too big to be stored without compression */
			/* Iterate over the data set to see how much memory we need to allocate */
			const size_t num_nks = mo.get_global_prof_by_name(gdb_prof_str_NK_DOMAIN).getiData().size();
			const size_t num_ks = mo.get_global_prof_by_name(gdb_prof_str_K_DOMAIN_CENTERED).getData().size();
			const int numfields = mo.get_global_param_by_name(gdb_str_NUMFIELDS).getiValue();	
			const rank4 bounds( num_nks, num_ks, numfields, numfields );
			size_t neededSize = 0;
			for(int t = 0; t <= maxTime; t++){
				for(int tprime = equalTimeOnly ? t : 0; tprime <= t; tprime++){
					const regset data = dataSet.getvar(t,tprime); // regsets are special pointers, actually create the object
					for( rank4iter rIter( bounds ); rIter < bounds; ++rIter ) {
						data.setprimary();
						const qgcomplex & primary_data = data[*rIter];
						data.setintermediate();
						const qgcomplex & intermediate_data = data[*rIter];
						if( shouldWriteComplexData(primary_data) || shouldWriteComplexData(intermediate_data) ){
							neededSize ++;	
						}

					}
				}	

			}
			printf("Needed size = %i\n", static_cast<int>(neededSize));
			complex<double> czero(0,0);	
			mxComplexRef initializer( czero );
			mxMatrix<uint8_T> mxIndices( neededSize, 6, -1 );
			mxVector< mxComplexRef > mxDataPrimary( neededSize, initializer );
			mxVector< mxComplexRef > mxDataIntermediate( neededSize, initializer );
			size_t ii = 0;
			for(int t = 0; t <= maxTime; t++){
				for(int tprime = equalTimeOnly ? t : 0; tprime <= t; tprime++){
					const regset data = dataSet.getvar(t,tprime); // regsets are special pointers, actually create the object
					for( rank4iter rIter( bounds ); rIter < bounds; ++rIter ) {
						data.setprimary();
						const qgcomplex & primary_data = data[*rIter];
						data.setintermediate();
						const qgcomplex & intermediate_data = data[*rIter];
						if( shouldWriteComplexData(primary_data) || shouldWriteComplexData(intermediate_data) ){
							if( ii >= neededSize ){
								printf("ii=%i exceeded neededSize=%i\n", static_cast<int>(ii), static_cast<int>(neededSize));
								continue;
							}
							mxIndices(ii,0) = t;
							mxIndices(ii,1) = tprime;
							for(size_t iter = 0; iter < 4; iter++){
								mxIndices(ii, iter+2) = static_cast<uint8_T>( (*rIter)[iter] ) + 1U;
							}
							mxDataPrimary[ii] = primary_data;
							mxDataIntermediate[ii] = intermediate_data;
							ii++;
						}
					}
				}
			}
			mxStruct mxSsystemField;
			mxSsystemField.setField( LABEL_MODEL_SSYSTEM_REGSET_INDICES, mxIndices );
			mxSsystemField.setField( LABEL_MODEL_SSYSTEM_REGSET_PRIMARY_DATA, mxDataPrimary );
			mxSsystemField.setField( LABEL_MODEL_SSYSTEM_REGSET_INTERMEDIATE_DATA, mxDataIntermediate );
			mxDest.setField( fieldName, mxSsystemField );
		}
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		void writeModelStateToFile() {
			mout.mxWriteVariable(modelState, LABEL_MODEL_STATE);
		}
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
		void execute() {
			checkMatFileStatus();
			saveModelParameters();
			saveModelProfiles();
			saveModelScriptInfo();
			saveModelDimensions();
			saveModelPSystem();
			saveModelSSystem();	
			writeModelStateToFile();
		}
		// - - - - - - - - - - - - - - - - - - - - - - - - - -
	public:
		mxSaveModelState(const model_class & model, const char * fileName) :
			mxModelInterface(fileName),
			mo(model),
			mout(fileName)
	{
		try {
			execute();
			success = true;
		}catch(std::exception & ex){
			dumpErrors();
			errorCleanUp();
		}
	};
};

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
// static functions that are called from outside
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
bool mx_save_model_state(const model_class & mo, const char * fname){
	return mxSaveModelState( mo, fname ) == true;
}
#endif
